################################################
#Code to reproduce Figure 2 
################################################

library(rio)
library(tidyverse)
library(ggstance)

#####################
#Functions
####################

`%notin%` <- Negate(`%in%`)

m1_func <- function(x) {
  x %>%
    filter(parm %notin% c("_cons", "1.post_soviet", "demyrs_01", "ethnic_01")) %>%
    filter(!is.na(z)) %>%
    mutate(model = "Model 1")
}

m2_func <- function(x) {
  x %>%
    filter(parm %notin% c("_cons", "1.post_soviet", "demyrs_01", "ethnic_01")) %>%
    filter(!is.na(z)) %>%
    mutate(model = "Model 2")
}

m3_func <- function(x) {
  x %>%
    filter(parm %notin% c("_cons", "1.post_soviet", "demyrs_01", "ethnic_01")) %>%
    filter(!is.na(z)) %>%
    mutate(model = "Model 3")
}

#####################
#Vote
####################

v1 <- import("vote_elections1.dta") %>% m1_func()
v2 <- import("vote_elections2.dta") %>% m2_func()
v3 <- import("vote_elections3.dta") %>% m3_func()

#####################
#Support Others
####################
s1 <- import("support_others1.dta") %>% m1_func()
s2 <- import("support_others2.dta") %>% m2_func()
s3 <- import("support_others3.dta") %>% m3_func()

#####################
#Obey
####################
o1 <- import("obey1.dta") %>% m1_func()
o2 <- import("obey2.dta") %>% m2_func()
o3 <- import("obey3.dta") %>% m3_func()

#####################
#Opinions
####################
op1 <- import("opinions1.dta") %>% m1_func()
op2 <- import("opinions2.dta") %>% m2_func()
op3 <- import("opinions3.dta") %>% m3_func()

#####################
#Volunteer
####################

vo1 <- import("volunteer1.dta") %>% m1_func()
vo2 <- import("volunteer2.dta") %>% m2_func()
vo3 <- import("volunteer3.dta") %>% m3_func()

#####################
#Be Active in Politics
####################

a1 <- import("beactive1.dta") %>% m1_func()
a2 <- import("beactive2.dta") %>% m2_func()
a3 <- import("beactive3.dta") %>% m3_func()

#####################
#Combine (Export)
####################

comb <- bind_rows(v1, v2, v3, 
                  a1, a2, a3, 
                  o1, o2, o3, 
                  op1, op2, op3, 
                  s1, s2, s3,
                  vo1, vo2, vo3)


#####################
#Figure
####################

comb1a <- comb %>%
  filter(model == "Model 3") %>%
  filter(parm %notin% c("2.edulvla", "3.edulvla", "4.edulvla", 
                        "0.edulvlfa", "2.edulvlfa", "3.edulvlfa", "4.edulvlfa", 
                        "3.born_when1")) %>%
  mutate(parm = factor(parm, 
                       levels=c("self_te", "cons_open", 
                                "age", "1.gender", 
                                "1.minority1", "1.union1", 
                                "5.edulvla", "5.edulvlfa", 
                                "2.born_when1", "socialtrust", "sclmeet", 
                                "religiosity", "yrlvdae"), 
                       labels=c("Self Transcendence", "Conservation", 
                                "Age", "Female", 
                                "Minority", "Union Mmbr.", 
                                "Tertiary Ed (Self)", "Tertiary Ed (Father)", 
                                "Came to Cntry <20yrs Ago", "Interpersonal Trust", 
                                "Social Engagement", "Religiosity", "Res. Stability")), 
         eq = factor(eq, 
                     levels=c("vote_elections", "beactive", "volunteer", 
                              "obey", "opinions", "support_others"), 
                     labels=c("Vote", "Be Active", "Volunteer", 
                              "Obey", "Opinions", "Support Others")), 
         values = ifelse(parm %in% c("Self Transcendence", "Conservation"), 1, 0), 
         values = factor(values, levels=c(1,0), 
                         labels=c("Values", "Controls")))
                       
                       

ggplot(comb1a, aes(x=estimate, y=parm)) + 
  geom_pointrangeh(aes(xmin=min95, xmax=max95)) + 
  facet_grid(values ~ eq, labeller = label_wrap_gen(width = 25), scales="free", space="free") + 
  geom_vline(xintercept=0, linetype="dashed", color="red") + 
  labs(y=NULL, x = "Coefficient") + 
  geom_text(aes(label=round(estimate, 2)), 
            nudge_y=0.20) +
  theme_bw(14) + 
  theme(strip.text.y = element_text(angle=0), 
        axis.title.y = element_text(margin = margin(r=15)), 
        strip.text = element_text(face = "bold")) 

ggsave("figure2.png", height=8, width=12, dpi=600)

                      